######
# This code reproduces the analysis in
# Ahlquist & Levi 2013 Ch. 7 (with Amanda Clayton) 
# LAST UPDATED: 19 August 2013
# UPDATED BY: ABC
# NOTE: this uses Amelia 1.6.3, ZeligChoice 0.7-0 and Zelig 4.1-2
######

rm(list=ls())
##load libraries 
library(MatchIt)
library(car)
library(vcd)
library(VGAM)
library(MASS)
library(Amelia)
library(foreign)
library(modeest)
library(ineq)


setwd("/Users/amandaclayton/Desktop/DataVerse")
modeldata<-read.csv("ALBookCh7.csv")
modeldata<-modeldata[,2:24]
names(modeldata)

##creating model data

##Use Amelia to impute missing data
bds1 <- matrix(c(4, 18, 96), nrow = 1, ncol = 3)  #set bounds for age (actual range in data)
bds2 <- matrix(c(2, 0, 52), nrow = 1, ncol = 3)  #set bounds for ilwu tenure (actual range in data)
bds3<-rbind(bds1, bds2)
 
m<-20
a.out <- amelia(modeldata, m = m,  ord=c("income", "edu", "round","nafta",
 "ilwu_fam"), nom=c("mode", "local", "imports_cat", "pol_party","white", "female","A", "B", "ilwu", "union_fam", "hispanic", "union", "vote", "protest", "contribute"), bounds=bds3)

##average across all 20 datasets
## no missingness in ILWU and Local, date
## mode for white, female, union family member and hispanic
## mean for age, education, income, nafta, ilwu tenure  

#create age average

temp<-0
for(i in 1:m){
  temp<-temp+a.out$imputations[[i]]$age}
comb.age<-temp/m

#create ilwu tenure average
temp<-0
for(i in 1:m){
  temp<-temp+a.out$imputations[[i]]$ilwu_years}
comb.ilwu_yrs<-round(temp/m,1)

#create education average
temp<-0
for(i in 1:m){
  temp<-temp+a.out$imputations[[i]]$edu}
comb.edu<-round(temp/m,0)

#create income average 
temp<-0
for(i in 1:m){
  temp<-temp+a.out$imputations[[i]]$income}
comb.income<-round(temp/m,0)

#create support for NAFTA average 
temp<-0
for(i in 1:m){
  temp<-temp+a.out$imputations[[i]]$nafta}
comb.nafta<-round(temp/m,0)

#pol party mode (used for specification checks)
temp<-NULL
for(i in 1:m){
  temp<-cbind(temp,a.out$imputations[[i]]$pol_party)}
party.mode<-apply(temp, 1, mfv)

#randomly filling in value where there are multiple modes
for(i in 1:length(party.mode)){
  if(sum(party.mode[[i]]==c(1,2))==2) party.mode[[i]]<-sample(c(1,2),1)
  if(sum(party.mode[[i]]==c(2,1))==2) party.mode[[i]]<-sample(c(1,2),1)
  if(sum(party.mode[[i]]==c(2,3))==2) party.mode[[i]]<-sample(c(3,2),1)
  if(sum(party.mode[[i]]==c(3,2))==2) party.mode[[i]]<-sample(c(3,2),1)
  if(sum(party.mode[[i]]==c(3,1))==2) party.mode[[i]]<-sample(c(3,1),1)
  if(sum(party.mode[[i]]==c(1,3))==2) party.mode[[i]]<-sample(c(3,1),1)
}
party.mode<-as.numeric(party.mode)


##union family member
temp<-NULL
for(i in 1:m){
  temp<-cbind(temp,a.out$imputations[[i]]$union_fam)}
union_fam.mode<-apply(temp, 1, mfv)

#ranomly filling in value where there are multiple modes
for(i in 1:length(union_fam.mode)){
  if(sum(union_fam.mode[[i]]==c(0,1))==2) union_fam.mode[[i]]<-sample(c(0,1),1)
  if(sum(union_fam.mode[[i]]==c(1,0))==2) union_fam.mode[[i]]<-sample(c(0,1),1)
}
union_fam.mode<-as.numeric(union_fam.mode)

##white
temp<-NULL
for(i in 1:m){
  temp<-cbind(temp,a.out$imputations[[i]]$white)}
white.mode<-apply(temp, 1, mfv)

#ranomly filling in value where there are multiple modes
for(i in 1:length(white.mode)){
  if(sum(white.mode[[i]]==c(0,1))==2) white.mode[[i]]<-sample(c(0,1),1)
  if(sum(white.mode[[i]]==c(1,0))==2) white.mode[[i]]<-sample(c(0,1),1)
}
white.mode<-as.numeric(white.mode)

###female
temp<-NULL
for(i in 1:m){
  temp<-cbind(temp,a.out$imputations[[i]]$female)}
female.mode<-apply(temp, 1, mfv)

#ranomly filling in value where there are multiple modes
for(i in 1:length(female.mode)){
  if(sum(female.mode[[i]]==c(0,1))==2) female.mode[[i]]<-sample(c(0,1),1)
  if(sum(female.mode[[i]]==c(1,0))==2) female.mode[[i]]<-sample(c(0,1),1)
}
female.mode<-as.numeric(female.mode)


##union member
temp<-NULL
for(i in 1:m){
  temp<-cbind(temp,a.out$imputations[[i]]$union)}
union.mode<-apply(temp, 1, mfv)

#randomly filling in value where there are multiple modes
for(i in 1:length(union.mode)){
  if(sum(union.mode[[i]]==c(0,1))==2) union.mode[[i]]<-sample(c(0,1),1)
  if(sum(union.mode[[i]]==c(1,0))==2) union.mode[[i]]<-sample(c(0,1),1)
}
union.mode<-as.numeric(union.mode)

##DV imports cat
temp<-NULL
for(i in 1:m){
  temp<-cbind(temp,a.out$imputations[[i]]$imports_cat)}
imports.mode<-apply(temp, 1, mfv)

for(i in 1:length(imports.mode)){
  if(sum(imports.mode[[i]]==c(1,2))==2) imports.mode[[i]]<-sample(c(1,2),1)
  if(sum(imports.mode[[i]]==c(2,1))==2) imports.mode[[i]]<-sample(c(1,2),1)
  if(sum(imports.mode[[i]]==c(2,3))==2) imports.mode[[i]]<-sample(c(3,2),1)
  if(sum(imports.mode[[i]]==c(3,2))==2) imports.mode[[i]]<-sample(c(3,2),1)
  if(sum(imports.mode[[i]]==c(3,1))==2) imports.mode[[i]]<-sample(c(3,1),1)
  if(sum(imports.mode[[i]]==c(1,3))==2) imports.mode[[i]]<-sample(c(3,1),1)
  if(sum(imports.mode[[i]]==c(0,1))==2) imports.mode[[i]]<-sample(c(0,1),1)
  if(sum(imports.mode[[i]]==c(1,0))==2) imports.mode[[i]]<-sample(c(0,1),1)
  if(sum(imports.mode[[i]]==c(0,2))==2) imports.mode[[i]]<-sample(c(0,2),1)
  if(sum(imports.mode[[i]]==c(2,0))==2) imports.mode[[i]]<-sample(c(0,2),1)
  if(sum(imports.mode[[i]]==c(0,3))==2) imports.mode[[i]]<-sample(c(0,3),1)
  if(sum(imports.mode[[i]]==c(3,0))==2) imports.mode[[i]]<-sample(c(0,3),1)
}

imports.mode<-as.numeric(imports.mode)

###hispanic
temp<-NULL
for(i in 1:m){
  temp<-cbind(temp,a.out$imputations[[i]]$hispanic)}
hispanic.mode<-apply(temp, 1, mfv)

#ranomly filling in value where there are multiple modes
for(i in 1:length(hispanic.mode)){
  if(sum(hispanic.mode[[i]]==c(0,1))==2) hispanic.mode[[i]]<-sample(c(0,1),1)
  if(sum(hispanic.mode[[i]]==c(1,0))==2) hispanic.mode[[i]]<-sample(c(0,1),1)
}
hispanic.mode<-as.numeric(hispanic.mode)

###A-status
temp<-NULL
for(i in 1:m){
  temp<-cbind(temp,a.out$imputations[[i]]$A)}
A.mode<-apply(temp, 1, mfv)

#ranomly filling in value where there are multiple modes
for(i in 1:length(A.mode)){
  if(sum(A.mode[[i]]==c(0,1))==2) A.mode[[i]]<-sample(c(0,1),1)
  if(sum(A.mode[[i]]==c(1,0))==2) A.mode[[i]]<-sample(c(0,1),1)
}
A.mode<-as.numeric(A.mode)

###B-status
temp<-NULL
for(i in 1:m){
  temp<-cbind(temp,a.out$imputations[[i]]$B)}
B.mode<-apply(temp, 1, mfv)

#randomly filling in value where there are multiple modes
for(i in 1:length(B.mode)){
  if(sum(B.mode[[i]]==c(0,1))==2) B.mode[[i]]<-sample(c(0,1),1)
  if(sum(B.mode[[i]]==c(1,0))==2) B.mode[[i]]<-sample(c(0,1),1)
}
B.mode<-as.numeric(B.mode)

##Protest
temp<-NULL
for(i in 1:m){
  temp<-cbind(temp,a.out$imputations[[i]]$protest)}
protest.mode<-apply(temp, 1, mfv)
for(i in 1:length(protest.mode)){
  if(sum(protest.mode[[i]]==c(0,1))==2) protest.mode[[i]]<-sample(c(0,1),1)
  if(sum(protest.mode[[i]]==c(1,0))==2) protest.mode[[i]]<-sample(c(0,1),1)
}

protest.mode<-as.numeric(protest.mode)

##Contributions
temp<-NULL
for(i in 1:m){
  temp<-cbind(temp,a.out$imputations[[i]]$contribute)}
contribute.mode<-apply(temp, 1, mfv)
for(i in 1:length(contribute.mode)){
  if(sum(contribute.mode[[i]]==c(0,1))==2) contribute.mode[[i]]<-sample(c(0,1),1)
  if(sum(contribute.mode[[i]]==c(1,0))==2) contribute.mode[[i]]<-sample(c(0,1),1)
}
contribute.mode<-as.numeric(contribute.mode)


##Vote 
temp<-NULL
for(i in 1:m){
  temp<-cbind(temp,a.out$imputations[[i]]$vote)}
vote.mode<-apply(temp, 1, mfv)
for(i in 1:length(vote.mode)){
  if(sum(vote.mode[[i]]==c(0,1))==2) vote.mode[[i]]<-sample(c(0,1),1)
  if(sum(vote.mode[[i]]==c(1,0))==2) vote.mode[[i]]<-sample(c(0,1),1)
}
vote.mode<-as.numeric(vote.mode)

comb.data<-as.data.frame(cbind(modeldata4$ilwu, modeldata4$local, round(a.out$imputations[[1]]$date, 0),
                               round(comb.age, 0), white.mode, comb.edu, female.mode, hispanic.mode, A.mode, B.mode, party.mode,
                               comb.income, a.out$imputations[[1]]$mode, union.mode, comb.ilwu_yrs, union_fam.mode, vote.mode,
                               contribute.mode, protest.mode,imports.mode , comb.nafta))

names(comb.data)<-c("ilwu", "local", "date", "age", "white", "edu", "female", "hispanic", "A", "B", "pol_party", "income", "mode",
                    "union", "ilwu_yrs", "union_fam", "vote", "contribute", "protest", "imports_cat","nafta" )

#recode factor variables as factors
comb.data$imports_cat<-as.factor(comb.data$imports_cat)
comb.data$pol_party<-as.factor(comb.data$pol_party)
comb.data$local<-as.factor(comb.data$local)
comb.data$mode<-as.factor(comb.data$mode)

###Matching
m.out1<-matchit(ilwu~white+as.factor(local)+age+edu+female+ date+hispanic,
                method="genetic", discard="both",data = comb.data, ratio=2)
match.out<-match.data(m.out1)
##matching diagnostics
summary(m.out1)
#plot(m.out1) 

#This code produces the plots in matching appendix
#This code produces Figure 7.12
plot(m.out1, type="hist") 

#This code produces Figure 7.13
plot(m.out1, type="jitter", interactive=F) 

#This code produces Figure 7.11
##100 x the mean diff/pooled SD 
summ.plot <- summary(m.out1, standardize=TRUE) 
plot(summ.plot, interactive=F, lines=NULL)
text(x=rep(0.9, 9), y=c(1.95, 1.55, 0.86, 0.74, 0.68, 0.35, 0.24, 0.158, 0.115, 0.065),
     labels=c("distance", "female", "age", "education", "Seattle", "Tacoma", "Long Beach",
       "Hispanic", "white", "date"), cex=0.7)


############# Trade opinion models for Ahlquist & Levi ch. 7 ###########

####This code produces Figure 7.2
##split data ILWU/RDD
ilwudata<-subset(modeldata4, ilwu==1) #638
rdddata<-subset(modeldata4, ilwu==0) #604
ilwu_unmatched<-table(ilwudata$imports_cat2)/length(na.omit(ilwudata$imports_cat2))
 #not imputed 
rdd_unmatched<-table(rdddata$imports_cat2)/length(na.omit(rdddata$imports_cat2))
 #not imputed 
##to compare w/ ANES data (from the 2004 ANES time series study)

nes<-read.dta("anes2004TS.dta")
imports<-nes$V045114
imports<-as.character(imports)
imports[imports=="1. Favor"]<-1 
imports[imports=="5. Oppose"]<-2 
imports[imports=="7. Haven't thought much about this"]<-3 
imports[imports=="8. Don't know"]<-NA 
imports[imports=="9. Refused"]<-NA
ANES<-table(imports)/(length(na.omit(imports)))

unmatched_table<-(rbind(ilwu_unmatched, rdd_unmatched, ANES))

#for matched data
##split matched data ILWU/RDD
MatchILWUdata<-subset(match.out, ilwu==1) #637
MatchRDDdata<-subset(match.out, ilwu==0) #256
ilwu_matched<-table(MatchILWUdata$imports_cat)/length(na.omit(MatchILWUdata$imports_cat))
rdd_matched<-table(MatchRDDdata$imports_cat)/length(na.omit(MatchRDDdata$imports_cat))

matched_table<-rbind(ilwu_matched, rdd_matched)

par(mfrow=c(1,2))
barplot(unmatched_table, beside=T, main="Support for Import Restrictions \n Unmatched Sample", ylim=c(0, 0.6),
        xlab=c(""), xaxt="n", col=c("black", "dark grey", "white"), width=c(0.2, 0.2), cex.main=0.7, space=c(0,1))
axis(side=1, at=c(0.35, 1.2, 2.2), labels=c("Favor", "Oppose", "Haven't Thought \n About It"),
     tick=FALSE, cex.axis=0.7)
legend("topleft", fill= c("black", "dark grey", "white"), border="black", col=c("black", "dark grey", "white"),
       legend=c("ILWU", "RDD", "ANES"), cex=0.6)

barplot(matched_table, beside=T, main="Support for Import Restrictions \n Matched Sample", ylim=c(0, 0.6), xlab=c(""),
        xaxt="n", col=c("black", "dark grey"), width=c(0.2, 0.2), cex.main=0.7, space=c(0,1))
axis(side=1, at=c(0.35, 1, 1.6), labels=c("Favor", "Oppose", "Haven't Thought \n About It"), tick=FALSE, cex.axis=0.7)
legend("topleft", fill=c("black", "dark grey"), border="black", col=c("black", "dark grey"),
       legend=c("ILWU", "RDD"), cex=0.6)

###This code produces Table 7.1
##main imports model
##NOTE: income, PID, survey mode all excluded as post-treatment variables
#
impres.1<-vglm(as.factor(imports_cat)~ ilwu + A + B +
              white + female + as.factor(local) +
              age + I(age^2) + edu + date + hispanic,
      family = multinomial,data=match.out, weights = match.out$weights) 
summary(impres.1)
yhats<-predictvglm(impres.1, type="response")
#with alternate reference category:
alt<-as.data.frame(match.out)
##original coding: 1=favor, 2=oppose, 3=dk (ref cat)
alt$imports_cat<-recode(alt$imports_cat, "1=2; 2=3; 3=1") #1=don't know, 2=favor, 3=oppose (ref cat)
impres.alt<-vglm(as.factor(imports_cat)~ ilwu + A + B +
              white + female + as.factor(local) +
              age + I(age^2) + edu + date + hispanic,
      family = multinomial,data=alt, weights = alt$weights) 
summary(impres.alt)

##Simulations and ternary plots
#calculating risk ratio 
x.trt.A<- c(1,1,0,1,0,3,mean(match.out$age), mean(match.out$edu),
            round(mean(match.out$date),0),0)
x.cntrl<- c(0,0,0,1,0,3,mean(match.out$age), mean(match.out$edu),
            round(mean(match.out$date),0),0)
x.trt.A<- data.frame(t(x.trt.A))
x.cntrl<- data.frame(t(x.cntrl))

names(x.trt.A)<-names(x.cntrl)<-c("ilwu", "A", "B", "white", "female", "local",
              "age", "edu" , "date" , "hispanic")
risk.rat<-predictvglm(impres.1, type="response", newdata=x.trt.A)/predictvglm(impres.1, type="response", newdata=x.cntrl)


y.hat.rdd<-yhats[match.data(m.out1)$ilwu==0,]
y.hat.ilwu<-yhats[match.data(m.out1)$ilwu==1,]
y.hat.A<-yhats[match.data(m.out1)$A==1,]
y.hat.B<-yhats[match.data(m.out1)$B==1,]
y.hat.C<-yhats[match.data(m.out1)$A==0 & match.data(m.out1)$B==0 & match.data(m.out1)$ilwu==1,]
y.hats.ABC<-rbind(y.hat.A,y.hat.B, y.hat.C)
A.men<-c(rep(1,nrow(y.hat.A)), rep(0,nrow(y.hat.B)), rep(0,nrow(y.hat.C)))
B.men<-c(rep(0,nrow(y.hat.A)), rep(1,nrow(y.hat.B)), rep(0,nrow(y.hat.C)))
y.hats.ABC<-cbind(y.hats.ABC,A.men, B.men)

##This code produces Figure 3
ternaryplot(
  yhats,
  pch = ifelse(match.out$ilwu == 0, 16, 2),
  col = ifelse(match.out$ilwu == 0, "black" , grey(0.5)),
  main="Predicted Probabilities \n ILWU/RDD Attitudes on Placing Limits on Imports",
  dimnames=c("favor","oppose","haven't thought about it"),
  cex = .7
)
grid_legend(0.2,0.63, pch=c(2,16), col=c(grey(0.5), "black"), labels=c("ILWU observations", "RDD observations"),
            title="", vgap=1, frame=F)

##This code produces Figure 7.4
###A men vs casuals

###A men vs casuals
ternaryplot(
   y.hats.ABC[,1:3],
  pch = ifelse(A.men==1 & B.men!=1, 17,
  		ifelse(A.men!=1 & B.men==1, 3, 
  		1)),
  col = ifelse(A.men==1 & B.men!=1, grey(0.7),
  		ifelse(A.men!=1 & B.men==1, grey(0.5), 
  		"black")),
  main="Predicted Probabilities \n A/B/Casuals' Attitudes on Placing Limits on Imports",
  dimnames=c("favor","oppose","haven't thought about it"),
  cex = .6
)


grid_legend(0.2,0.63, pch=c(17,3,1), col=c(grey(0.7), grey(0.5), "black"), labels=c("Class-A", "Class-B", "Casuals"),
            title="", vgap=1, frame=F)

#### compare NAFTA responses to WSV data
#this code produces Figure 7.6
setwd("/Users/amandaclayton/Documents/labor/trade attitudes/other surveys")
wvs<-read.dta("wvs.dta") ##from 2006, most recent tabulated year of WVS
wvsUS<-as.data.frame(subset(wvs, wvs$s024a=="united states (5)"))

nafta<-wvsUS$e069_24
table(nafta)
nafta1<-recode(nafta, "'a great deal' = 1; 'quite a lot' = 2; 'not very much' = 3; 'none at all' = 4; NA=NA")
wvs_nafta<-table(nafta1)/length(na.omit(nafta1))  #n = 1170
ilwu_nafta<-table(ilwudata$nafta)/length(na.omit(ilwudata$nafta))
 #not imputed 
rdd_nafta<-table(rdddata$nafta)/length(na.omit(rdddata$nafta))
 #not imputed 

naftatable<-rbind(ilwu_nafta, rdd_nafta, wvs_nafta)

par(mfrow=c(1,1))

barplot(naftatable, beside=T, main="Respondents' 'Confidence' in NAFTA",ylim=c(0, 0.7),
        xlab=c(""), xaxt="n", col=c("black", "dark grey", "white"), width=c(0.2, 0.2), cex.main=1, space=c(0,1))
axis(side=1, at=c(0.5, 1.3, 2.1, 2.9), labels=c("A great deal", "Quite a lot", "Not very much", "None at all"),
     tick=FALSE, cex.axis=0.8)
legend("topleft", #pch=15, 
fill=c("black", "dark grey", "white"), border="black", col=c("black", "dark grey", "white"),
       legend=c("ILWU", "RDD", "WVS"), cex=0.8)


##This code produces Table 7.2
##main NAFTA model
library(Zelig)
library(ZeligChoice)
nafta.fit<-zelig(as.factor(nafta)~ilwu + A + B + white + female +
              as.factor(local) + age + I(age^2) + edu + date + hispanic, data=match.data(m.out1),
              weights = match.data(m.out1)$weights, model="oprobit") 
summary(nafta.fit)

##1=a great deal, 2=quite a lot, 3=not very much, 4=none at all

#Simulations and first difference plots
#This code produces Figure 7.7 
##RDD vs ILWU 
x.ILWU <- setx(nafta.fit, ilwu=1, A=1, B=0, white=1, female = 0, hispanic=0)  
x.RDD <- setx(nafta.fit, ilwu = 0, A=0, B=0, white=1, female = 0, hispanic=0)
x.cas <- setx(nafta.fit, ilwu = 1, A=0, B=0, white=1, female = 0, hispanic=0)

s.out.1<-sim(nafta.fit, x=x.ILWU, x1=x.RDD) 
s.out.2<-sim(nafta.fit, x=x.ILWU, x1=x.cas) 

plot(density(s.out.1$qi$fd[,4], adjust = 2), xlab="probability", xlim=c(-.1,.45),
     ylim=c(0,15), lwd=1.5,
     main = "Implied effect on NAFTA confidence, ordered probit")
lines(density(s.out.2$qi$fd[,4], adjust = 2), col="blue", lwd=1.5, lty=2)

text(0.21,14.8,"Pr(none | Class-A) - Pr(none | Casual)", cex=.7)
text(.36,11,"Pr(none | Class-A) - Pr(none | RDD)", cex=.7)
abline(v=0, lty=3, col=grey(.3))



######
## matching by age and union tenure; comparing A/B attitutudes toward trade 
## This code produces Figure 7.5
AB<-subset(comb.data,A==1 | B==1)

m.out2<-matchit(A~age + female + white + hispanic + edu + mode, method="genetic", discard="both", data=AB)
m.out3<-matchit(A~ilwu_yrs + female + white + hispanic + edu + mode, method="genetic", discard="both", data=AB)
summary(m.out2)
summary(m.out3)


#difference by age
summary(match.data(m.out2))
A2<-subset(match.data(m.out2), A==1)
B2<-subset(match.data(m.out2), B==1)

age_match<-rbind(table(A2$imports_cat)/(nrow(A2)), table(B2$imports_cat)/(nrow(B2)))

par(mfrow=c(1,2))

barplot(age_match, beside=T, main="Support for Import Restrictions \n  Class-A/Class-B Members Matched on Age",
        ylim=c(0, 0.6), xlab=c(""), ylab=c("as % of matched Class-A and Class-B ILWU respondents"), xaxt="n",
        col=c("black", "dark grey"), width=c(0.2, 0.2), cex.main=0.7, space=c(0,1)) #, cex.axis=0.6)
#axis(2,cex.axis=0.2)
axis(side=1, at=c(0.35, 1, 1.7), labels=c("Favor", "Oppose", "Don't Know"), tick=FALSE, cex.axis=0.7)
legend("topleft", pch=15, col=c("black", "dark grey"), legend=c("A-member", "B-member"), cex=0.6)

#difference by union tenure
A3<-subset(match.data(m.out3), A==1)
B3<-subset(match.data(m.out3), B==1)

tenure_match<-rbind(table(A3$imports_cat)/(nrow(A3)),table(B3$imports_cat)/(nrow(B3)))

barplot(tenure_match, beside=T, main="Support for Import Restrictions \n Class-A/Class-B Members Matched on Union Tenure",
        ylim=c(0, 0.6), xlab=c(""), xaxt="n", col=c("black", "dark grey"), width=c(0.2, 0.2), cex.main=0.7, space=c(0,1))
axis(side=1, at=c(0.35, 1, 1.7), labels=c("Favor", "Oppose", "Don't Know"), tick=FALSE, cex.axis=0.7)
legend("topleft", pch=15, col=c("black", "dark grey"), legend=c("A-member", "B-member"), cex=0.6)



##Additional specifications for Online Appendix
###including partyID, income, & survey mode as covariates in the structural model
trade.fit.pidinc<-zelig(as.factor(imports_cat)~ilwu + ilwu:A + ilwu:B + white + female + as.factor(local) +
                        age + I(age^2) + edu + date + as.factor(mode) + hispanic + pol_party + income,
                        data=match.data(m.out1), weights = match.data(m.out1)$weights, model="mlogit") 

nafta.fit.pidinc<-zelig(as.factor(nafta)~ilwu + ilwu:A + ilwu:B + white + female + as.factor(local) +
                        age + I(age^2) + edu + date + as.factor(mode) + hispanic + pol_party + income,
                        data=match.data(m.out1), weights = match.data(m.out1)$weights, model="oprobit") 
summary(trade.fit.pidinc)
summary(nafta.fit.pidinc)
  
### Other matching exercises: matching on union family member and using union-only subsample

##### Matching on union family member
m.out5<-matchit(ilwu~white+as.factor(local)+age+edu+female+ date+hispanic+union_fam,
                method="genetic", discard="both",data=comb.data, ratio=2)
summary(m.out5)
###This code produces Table A2.1
##imports model
z.out5<-zelig(as.factor(imports_cat)~ilwu+ A + B +
              white+female+as.factor(local)+age+I(age^2)+edu+date+hispanic+union_fam,
              data=match.data(m.out5), weights = match.data(m.out5)$weights, model="mlogit") 
summary(z.out5)
#with alternate reference category:
alt5<-as.data.frame(match.data(m.out5))
##original coding: 1=favor, 2=oppose, 3=dk (ref cat)
alt5$imports_cat<-recode(alt5$imports_cat, "1=2; 2=3; 3=1") #1=don't know, 2=favor, 3=oppose (ref cat)

z.out6<-zelig(as.factor(imports_cat)~ilwu+ ilwu:A + ilwu:B +
              white+female+as.factor(local)+age+I(age^2)+edu+date+hispanic+union_fam,
              data=alt5, weights = alt5$weights, model="mlogit") 
summary(z.out6)

#This code produces Table A2.2
#nafta model
z.out7<-zelig(as.factor(nafta)~ilwu + ilwu:A + ilwu:B + white + female +
              as.factor(local) + age + I(age^2) + edu + date + hispanic + union_fam,
              data=match.data(m.out5), weights = match.data(m.out5)$weights, model="oprobit") 

summary(z.out7)
##1=a great deal, 2=quite a lot, 3=not very much, 4=none at all

#union-only subsample
union.sub<-subset(match.data(m.out1), union==1)

###This code produces Table A2.3
##imports model
z.out8<-zelig(as.factor(imports_cat) ~ ilwu + ilwu:A + ilwu:B + white + female + as.factor(local) +
              age + I(age^2) + edu + date+ hispanic,
              data=union.sub, weights = union.sub$weights,model="mlogit") 
summary(z.out8)

x1 <- setx(z.out8, ilwu=1, A=1, B=0, white=1, female = 0, hispanic = 0)  
xx <- setx(z.out8, ilwu=0, A=0, B=0, white=1, female = 0, hispanic = 0)

s.out.fd<-sim(z.out8, x=xx, x1=x1)
summary(s.out.fd) #getting risk ratios for paper

#with alternate reference category:
union.sub.alt<-union.sub
##original coding: 1=favor, 2=oppose, 3=dk (ref cat)
union.sub.alt$imports_cat<-recode(union.sub.alt$imports_cat, "1=2; 2=3; 3=1") #1=don't know, 2=favor, 3=oppose (ref cat)

z.out9<-zelig(as.factor(imports_cat) ~ ilwu+ ilwu:A + ilwu:B + white+female+as.factor(local)+
              age+I(age^2) + edu + date +hispanic,
              data=union.sub.alt, weights = union.sub.alt$weights, model="mlogit") 
summary(z.out9)

###This code produces Table A2.4
#nafta model
z.out10<-zelig(as.factor(nafta)~ ilwu+ ilwu:A + ilwu:B + white + female + as.factor(local) +
               age + I(age^2)+edu+date+ hispanic,
               data=union.sub.alt, weights = union.sub.alt$weights, model="oprobit") 
summary(z.out10)

##1=a great deal, 2=quite a lot, 3=not very much, 4=none at all

#####Political Behavior Models for Ahlquist & Levi ch. 7

####  PROTEST
unionmatch1<-subset(match.data(m.out1), ilwu==1) 
table1<-table(unionmatch1$protest)
rddmatch1<-subset(match.data(m.out1), ilwu==0) 
table2<-table(rddmatch1$protest)

protest_tab<-rbind((table1/nrow(unionmatch1)),(table2/nrow(rddmatch1)))

##This code produces Figure 7.8

par(mfrow=c(1,2))
barplot(protest_tab, beside=T, main="Total Respondents' Participation in a Political \n Demonstration in the Last 12 Months",
        ylim=c(0,1), xlab=c(""), ylab=c("as % of total respondents"), xaxt="n", col=c("black", "dark grey"), width=0.5, cex.main=0.75)
axis(side=1, at=c(1, 2.5), labels=c("No", "Yes"), tick=FALSE, cex.axis=1)
legend("topright", pch=15, col=c("black", "dark grey"), legend=c("ILWU", "RDD"), cex=0.6)

##by rank
A1<-subset(unionmatch1,A==1)
B1<-subset(unionmatch1,B==1)
C1<-subset(unionmatch1,A==0 & B==0)
tableA<-table(A1$protest)
tableB<-table(B1$protest)
tableC<-table(C1$protest)

protest_rank<-rbind((tableA/nrow(A1)),(tableB/nrow(B1)), (tableC/nrow(C1)))

barplot(protest_rank, beside=T, main="ILWU Respondents' Participation in a Political \n Demonstration in the Last 12 Months",
        ylim=c(0,1), xlab=c(""), ylab=c("as % of total ILWU respondents"), xaxt="n", col=c("black", "dark grey", "white"), width=0.5, cex.main=0.75)
axis(side=1, at=c(1.25, 3.25), labels=c("No", "Yes"), tick=FALSE, cex.axis=1)
legend("topright", #pch=15, 
fill=c("black", "dark grey", "white"), border="black", col=c("black", "dark grey", "white"), legend=c("A-member", "B-member", "Casual"), cex=0.6)

#Protest logistic regression model
## Ahlquist & Levi table 7.3
protest.fit<-zelig(as.factor(protest)~ ilwu+ A + B +white+hispanic+female+as.factor(local)+as.numeric(age)+ I(age^2) +
                   as.numeric(edu)+date, data=match.data(m.out1), weights = match.data(m.out1)$weights, model="logit") 
summary(protest.fit)
## interpretation
x.ILWU <- setx(protest.fit, ilwu=1, A=1, B=0, white=1, hispanic=0, female=0, edu=4)  
x.RDD <- setx(protest.fit, ilwu = 0, A=0, B=0, white=1, hispanic=0, female=0, edu=4)

s.protest<-sim(protest.fit, x=x.RDD, x1=x.ILWU)
summary(s.protest)

x.A <- setx(protest.fit, ilwu=1, A=1, B=0, white=1, hispanic=0, female=0,  edu=4) 
x.C <- setx(protest.fit, ilwu =1, A=0, B=0, white=1, hispanic=0, female=0,  edu=4)
s.protest<-sim(z.out1, x=x.C, x1=x.A)
summary(s.protest)

##### CONTRIBUTIONS
table1<-table(unionmatch1$contribute)
table2<-table(rddmatch1$contribute)

contribution_tab<-rbind((table1/nrow(unionmatch1)),(table2/nrow(rddmatch1)))

#This code produces Figure 7.9

par(mfrow=c(1,2))
barplot(contribution_tab, beside=T, main="Total Respondents' Contribution to a \n Political Candidate in the Last 12 Months",ylim=c(0,1), xlab=c(""), ylab=c("as % of total respondents"), xaxt="n", col=c("black", "dark grey"), width=0.5, cex.main=0.8)
axis(side=1, at=c(1, 2.5), labels=c("No", "Yes"), tick=FALSE, cex.axis=1)
legend("topright", pch=15, col=c("black", "dark grey"), legend=c("ILWU", "RDD"), cex=0.6)

##by rank
tableA<-table(A1$contribute)
tableB<-table(B1$contribute)
tableC<-table(C1$contribute)

contribution_rank<-rbind((tableA/nrow(A1)),(tableB/nrow(B1)), (tableC/nrow(C1)))

barplot(contribution_rank, beside=T, main="ILWU Respondents' Contribution to a \n Political Candidate in the Last 12 Months",
        ylim=c(0,1), xlab=c(""), ylab=c("as % of total ILWU respondents"), xaxt="n", col=c("black", "dark grey", "white"), width=0.5, cex.main=0.8)
axis(side=1, at=c(1.25, 3.25), labels=c("No", "Yes"), tick=FALSE, cex.axis=1)
legend("topright", #pch=15, 
fill=c("black", "dark grey", "white"), border="black", col=c("black", "dark grey", "white"), legend=c("A-member", "B-member", "Casual"), cex=0.6)


match.data<-as.data.frame(match.data(m.out1))
#baseline = referring to odd years - those that responded during 6/19/06 - 11/06/06; 12/1/07 - 11/3/08; 12/01/09 - 11/1/10
#referring to 06 midterms: 11/07/06 - 11/30/07 (187:370)
#referring to 08 prez: 11/4/08 - 11/30/09 (885:1026)
#referring to 10 midterms: 11/2/10 - 7/29/11 (1600:1867)

match.data$mid06<-recode(match.data$date, "187:370=1; else=0")
match.data$prez08<-recode(match.data$date, "885:1026=1; else=0")
match.data$mid10<-recode(match.data$date, "1600:1867=1; else=0")

#This code produces Table 7.4
contribute.fit<-zelig(as.factor(contribute)~ ilwu + A + B +white+hispanic+female+as.factor(local)+as.numeric(age)+ I(age^2)+as.numeric(edu)
              +date+mid06+prez08+mid10, data=match.data, weights = match.data$weights, model="logit") 
summary(contribute.fit)

##interpreting results 
##RDD vs ILWU 
x.ILWU <- setx(contribute.fit, ilwu=1, A=1, B=0, white=1, hispanic=0, female=0, edu=4, mid10=1, mid06=0, prez08=0) #income mean in matched data  
x.RDD <- setx(contribute.fit, ilwu = 0, A=0, B=0, white=1, hispanic=0, female=0, edu=4, mid10=1, mid06=0, prez08=0)

contribute.s<-sim(contribute.fit, x=x.RDD, x1=x.ILWU)
summary(contribute.s)

###A-men vs casuals
x.A <- setx(contribute.fit, ilwu=1, A=1, B=0, white=1, female=0, edu=4, mid10=1, mid06=0, prez08=0) 
x.C <- setx(contribute.fit, ilwu =1, A=0, B=0, white=1, female=0, edu=4, mid10=1, mid06=0, prez08=0)

contribute.s<-sim(contribute.fit, x=x.C, x1=x.A)
summary(contribute.s)

#####VOTE

table1<-table(unionmatch1$vote)
table2<-table(rddmatch1$vote)

vote_tab<-rbind((table1/nrow(unionmatch1)),(table2/nrow(rddmatch1)))

#This code produces Figure 7.10
par(mfrow=c(1,2))
barplot(vote_tab, beside=T, main="Whether Total Respondents Voted \n in the Last Election",ylim=c(0,1),
        xlab=c(""), ylab=c("as % of total respondents"), xaxt="n", col=c("black", "dark grey"), width=0.5, cex.main=0.8)
axis(side=1, at=c(1, 2.5), labels=c("No", "Yes"), tick=FALSE, cex.axis=1)
legend("topleft", pch=15, col=c("black", "dark grey"), legend=c("ILWU", "RDD"), cex=0.6)

##by rank
tableA<-table(A1$vote)
tableB<-table(B1$vote)
tableC<-table(C1$vote)

vote_rank<-rbind( (tableA/nrow(A1)),(tableB/nrow(B1)), (tableC/nrow(C1)))

barplot(vote_rank, beside=T, main="Whether ILWU Respondents Voted \n in the Last Election",ylim=c(0,1), xlab=c(""),
        ylab=c("as % of total ILWU respondents"), xaxt="n", col=c("black", "dark grey", "white"), width=0.5, cex.main=0.8)
axis(side=1, at=c(1.25, 3.25), labels=c("No", "Yes"), tick=FALSE, cex.axis=1)
legend("topleft", #pch=15, 
fill=c("black", "dark grey", "white"), border="black", col=c("black", "dark grey", "white"), legend=c("A-member", "B-member", "Casual"), cex=0.6)


#This code produces Table 7.5
vote.fit<-zelig(as.factor(vote)~ilwu+ A+ B +white+hispanic+female+as.factor(local)+as.numeric(age)+ I(age^2)+as.numeric(edu)+date+
                mid06+prez08+mid10, data=match.data, weights = match.data$weights, model="logit") 

summary(vote.fit)

##interpreting results 
##RDD vs ILWU 
x.ILWU <- setx(vote.fit, ilwu=1, A=1, B=0, white=1, hispanic=0,  female=0, edu=4, mid10=1, mid06=0, prez08=0)  
x.RDD <- setx(vote.fit, ilwu = 0, A=0, B=0, white=1, hispanic=0, female=0, edu=4, mid10=1, mid06=0, prez08=0)

vote.s<-sim(vote.fit, x=x.RDD, x1=x.ILWU) 
summary(vote.s) 

###A-men vs casuals
x.A <- setx(vote.fit, ilwu=1, A=1, B=0, white=1, hispanic=0, female=0, edu=4, mid10=1, mid06=0, prez08=0) 
x.C <- setx(vote.fit, ilwu =1, A=0, B=0, white=1, hispanic=0, female=0, edu=4, mid10=1, mid06=0, prez08=0)
vote.s<-sim(vote.fit, x=x.C, x1=x.A)
summary(vote.s)

### Additional matches and robustness for web appendix
#including having a union family member as a covariate
protest.fit.uf<-zelig(as.factor(protest)~ ilwu+ A + B + white+hispanic+female+as.factor(local)+as.numeric(age)+ I(age^2)+
                      as.numeric(edu)+date+union_fam, data=match.data(m.out5), weights = match.data(m.out5)$weights, model="logit") 
summary(protest.fit.uf)

match.data<-as.data.frame(match.data(m.out5))
match.data$mid06<-recode(match.data$date, "187:370=1; else=0")
match.data$prez08<-recode(match.data$date, "885:1026=1; else=0")
match.data$mid10<-recode(match.data$date, "1600:1867=1; else=0")

contribute.fit.uf<-zelig(as.factor(contribute)~ ilwu+ A + B + white+hispanic+female+as.factor(local)+as.numeric(age)+ I(age^2)+
                      as.numeric(edu)+date+union_fam +mid06+prez08+mid10, data=match.data, weights = match.data$weights, model="logit")
summary(contribute.fit.uf)

vote.fit.uf<-zelig(as.factor(vote)~ ilwu+ A + B + white+hispanic+female+as.factor(local)+as.numeric(age)+ I(age^2)+
                      as.numeric(edu)+date+union_fam +mid06+prez08+mid10, data=match.data, weights = match.data$weights, model="logit")
summary(vote.fit.uf)

#union-only subsample
protest.fit.us<-zelig(as.factor(protest)~ ilwu+ A + B + white+hispanic+female+as.factor(local)+as.numeric(age)+ I(age^2)+
                      as.numeric(edu)+date+union_fam, data=union.sub, weights = union.sub$weights, model="logit")
summary(protest.fit.us)

union.sub$mid06<-recode(union.sub$date, "187:370=1; else=0")
union.sub$prez08<-recode(union.sub$date, "885:1026=1; else=0")
union.sub$mid10<-recode(union.sub$date, "1600:1867=1; else=0")

contribute.fit.us<-zelig(as.factor(contribute)~ ilwu+ A + B + white+hispanic+female+as.factor(local)+as.numeric(age)+ I(age^2)+
                      as.numeric(edu)+date+union_fam +mid06+prez08+mid10, data=union.sub, weights = union.sub$weights, model="logit")
summary(contribute.fit.us)

vote.fit.us<-zelig(as.factor(vote)~ ilwu+ A + B + white+hispanic+female+as.factor(local)+as.numeric(age)+ I(age^2)+
                      as.numeric(edu)+date+union_fam +mid06+prez08+mid10, data=union.sub, weights = union.sub$weights, model="logit")
summary(vote.fit.us)
###END
